!!S Main
	program matrix
	
	implicit none
	
	integer i,j,ticks1,ticks2
	extended a(50,50), b(50,50), c(50,50)
	extended time1,time2
	
	integer ticks
	
	type *,'Matrix multiplication benchmark'
	type *,'-------------------------------'
	type *
	type *,'This program compares the number crunching power'
	type *,'of some of the popular MPW compilers.'
	type *,'Written under MPW 3.0 by J. Langowski / MacTutor 1989'
	type *
	type *,'Setting up 50x50 matrices...'

	ticks1 = ticks()
	
	do i=1,50
		do j=1,50
			a(i,j) = (i-1) + j-1
			b(j,i) = a(i,j)
		end do
	end do
	
	ticks2 = ticks()
	time1 = (ticks2-ticks1)/60.	
	type *
	write (6,'(f8.4,'' seconds for setting up matrices'')') time1
	
	ticks1 = ticks()
	call mat_mult_for3(c,50,a,50,b,50,50,50,50)
	
	ticks2 = ticks()
	time1 = (ticks2-ticks1)/60.	
	type *
	write (6,'(f8.4,'' seconds for multiplying matrices'',
     *	'' using FORTRAN routine, opt=3'')') time1
	type *,'c(25,25) = ',c(25,25)
	
	ticks1 = ticks()
	call mat_mult_for(c,50,a,50,b,50,50,50,50)
	
	ticks2 = ticks()
	time1 = (ticks2-ticks1)/60.	
	type *
	write (6,'(f8.4,'' seconds for multiplying matrices'',
     *	'' using FORTRAN routine, opt=0'')') time1
	type *,'c(25,25) = ',c(25,25)

	ticks1 = ticks()
	call mat_mult_for1(c,50,a,50,b,50,50,50,50)
	
	ticks2 = ticks()
	time1 = (ticks2-ticks1)/60.	
	type *
	write (6,'(f8.4,'' seconds for multiplying matrices'',
     *	'' using FORTRAN routine, hand-optimized'')') time1
	type *,'c(25,25) = ',c(25,25)
	
	ticks1 = ticks()
	call mat_mul_pas(c,%val(50),a,%val(50),b,%val(50),%val(50),%val(50),%val(50))
	
	ticks2 = ticks()
	time1 = (ticks2-ticks1)/60.	
	type *
	write (6,'(f8.4,'' seconds for multiplying matrices'',
     *	'' using PASCAL routine'')') time1
	type *,'c(25,25) = ',c(25,25)
	
	ticks1 = ticks()
	call mat_mul_pas_opt(c,%val(50),a,%val(50),b,%val(50),%val(50),%val(50),%val(50))
	
	ticks2 = ticks()
	time1 = (ticks2-ticks1)/60.	
	type *
	write (6,'(f8.4,'' seconds for multiplying matrices'',
     *	'' using PASCAL routine, hand-optimized'')') time1
	type *,'c(25,25) = ',c(25,25)
	
	ticks1 = ticks()
	call mat_mul_c(c,%val(50),a,%val(50),b,%val(50),%val(50),%val(50),%val(50))
	
	ticks2 = ticks()
	time1 = (ticks2-ticks1)/60.	
	type *
	write (6,'(f8.4,'' seconds for multiplying matrices'',
     *	'' using C routine'')') time1
	type *,'c(25,25) = ',c(25,25)
	
	ticks1 = ticks()
	call mat_mul_c_opt(c,%val(50),a,%val(50),b,%val(50),%val(50),%val(50),%val(50))
	
	ticks2 = ticks()
	time1 = (ticks2-ticks1)/60.	
	type *
	write (6,'(f8.4,'' seconds for multiplying matrices'',
     *	'' using C routine, hand-optimized'')') time1
	type *,'c(25,25) = ',c(25,25)
	
	end


!!S Main
	subroutine mat_mult_for3(c,nc,a,na,b,nb,n1,n2,n3)
c sets c=a.b
c na,nb,nc are first dimensions
c n1 n2 n3 are problem dimensions
c c is n1xn3
c a    n1 n2
c b    n2 n3
	implicit none
	
	integer na,nb,nc,n1,n2,n3
	integer*2 i,j,k
	extended c(nc,n3),a(na,n2),b(nb,n3)
	
	do k=1,n3
		do i=1,n1
			c(i,k) = 0x0
			do j=1,n2
				c(i,k) = c(i,k)+a(i,j)*b(j,k)
			end do
		end do
	end do
	return
	end

!!S Main
	subroutine mat_mult_for1(c,nc,a,na,b,nb,n1,n2,n3)
c sets c=a.b
c na,nb,nc are first dimensions
c n1 n2 n3 are problem dimensions
c c is n1xn3
c a    n1 n2
c b    n2 n3
	implicit none
	
	integer na,nb,nc,n1,n2,n3
	integer*2 i,j,k
	extended c(nc,n3),a(na,n2),b(nb,n3)
	extended sum 
	
	do k=1,n3
		do i=1,n1
			sum = 0x0
			do j=1,n2
				sum = sum+a(i,j)*b(j,k)
			end do
			c(i,k) = sum
		end do
	end do
	return
	end

!!S Main
	integer function ticks
	ticks = long(362)
	return
	end